#include "structMeshFunctions.hpp"
#include "macros.hpp"
#include <iostream>

namespace hms
{

template<typename MeshType, typename Callable1, typename Callable2, typename ... Ts>
void traverseEdgesSecondOrder(
	[[maybe_unused]] EdgeTraversal traversal,
	const MeshBase<MeshType>& mesh,
	const Array<bool, 4, 1>& meshPartAdjacency,
	[[maybe_unused]]Callable1 funcCentral,
	[[maybe_unused]]Callable2 funcMUSCL,
	Ts&&... args
){

	static_assert( std::is_invocable_r_v
		<void, Callable1, const MeshType&, Index, Ts ...>,
		"Callable1 must be invocable with const MeshType& as first "
		"and Index as second argument!"
	);
	static_assert( std::is_invocable_r_v
		<void, Callable2, const MeshType&, Index, Ts ...>,
		"Callable2 must be invocable with const MeshType& as first "
		"and Index as second argument!"
	);

	// Just a placeholder handling. This needs to be changed.
	if( (meshPartAdjacency == false).any() ){
		std::cout << "Please handle adjacency of meshparts.\n";
	}

	using T = EdgeTraversal;

	#define FORWARD_ARGS \
		mesh, funcCentral, funcMUSCL, std::forward<Ts>(args) ...

	if constexpr ( std::is_base_of_v<StructMeshBase<MeshType>, MeshType> ){
		switch (traversal){
			case T::ordered_secondOrder:
				orderedTraversal_secondOrder( FORWARD_ARGS );
				break;

			case T::rowwise_secondOrder:
				rowwiseTraversal_secondOrder( FORWARD_ARGS );
				break;

			case T::rowwise_secondOrder_seamed:
				rowwiseTraversal_secondOrder_seamed( FORWARD_ARGS );
				break;

			case T::blockwise_secondOrder:
				blockwiseTraversal_secondOrder( FORWARD_ARGS );
				break;

			case T::blockwise_secondOrder_seamed:
				blockwiseTraversal_secondOrder_seamed( FORWARD_ARGS );
				break;
			
			default:
				std::cerr << "WARNING: Unknown edge traversal type. "
					"Defaulting to ordered second order traversal.\n";
				orderedTraversal_secondOrder( FORWARD_ARGS );
				break;
		}
	}
	else {}
}

// traversal methods for higher order reconstruction schemes
TRAVERSAL_FUNC_SECOND_ORDER(orderedTraversal_secondOrder){
	STRUCTURED_TRAVERSAL_ASSERT;

	Index nx { mesh.derived().nCellsX() };
	Index ny { mesh.derived().nCellsY() };

	Index nh { nHorizontalEdges(mesh.derived().size()) };

	// INNER EDGES
	Index horiBeg { nx*2 };
	Index horiEnd { nh - horiBeg };

	Index nVertEdgesX { nx +1 };
	Index vertBeg { 2 }; 	// first vertical edge in a row that is neither 
							// boundary edge or direct neighbour to boundary 
							// edge of left boundary
							
	Index vertEnd { nVertEdgesX - 2}; 	// last vertical edge in a row that is neither 
									// boundary edge or direct neighbour to boundary 
									// edge of right boundary
	OMP_PRAGMA(parallel)
	{
		// horizontal edges - central
		OMP_PRAGMA(for schedule (static) nowait)
		for (Index e = 0; e < 2*nx; e++ ) {
			funcCentral( mesh, e, std::forward<Ts>(args) ... );
		}

		// horizontal edges - muscl
		OMP_PRAGMA(for schedule (static) nowait)
		for(Index e = horiBeg; e< horiEnd; e++){
			funcMUSCL( mesh, e, std::forward<Ts>(args) ... );
		}

		// last horizontal edges - central again
		OMP_PRAGMA(for schedule (static) nowait)
		for (Index e = (nh - 2*nx); e < nh; e++ ) {
			funcCentral( mesh, e, std::forward<Ts>(args) ... );
		}

		// loop for vertical edges.
		OMP_PRAGMA(for schedule (static) nowait)
		for( Index row = 0; row < ny; row++ ){ // loop through rows

			Index eFirst { nh + row * nVertEdgesX };
			Index eLast  { nh + row * nVertEdgesX + nx };
			
			// first two edges in a row are handled first order.
			funcCentral( mesh, eFirst, std::forward<Ts>(args) ... );
			funcCentral( mesh, eFirst+1, std::forward<Ts>(args) ... );

			// loop through edges in between the first two and last two in a row.
			for(Index i = vertBeg; i < vertEnd; i++) {
				Index e = nh + row * (nVertEdgesX) + i;
				funcMUSCL( mesh, e, std::forward<Ts>(args) ... ) ;
			}
			
			// last two edges in a row are handled first order again.
			funcCentral( mesh, eLast-1, std::forward<Ts>(args) ... );
			funcCentral( mesh, eLast, std::forward<Ts>(args) ... );
		}
	}
}

TRAVERSAL_FUNC_SECOND_ORDER(rowwiseTraversal_secondOrder){

	STRUCTURED_TRAVERSAL_ASSERT;
	auto callCentral = [&](Index e){
		funcCentral(mesh, e, std::forward<Ts>(args) ... );
	};
	auto callMUSCL = [&](Index e){
		funcMUSCL(mesh, e, std::forward<Ts>(args) ... );
	};

	Index
		nx { mesh.derived().size().x() },
		ny { mesh.derived().size().y() };
	
	Index nh { nHorizontalEdges( mesh.derived().size() ) };

	OMP_PRAGMA(parallel)
	{
		OMP_PRAGMA(for schedule(static) nowait)
		// loop through first two rows
		for ( int row = 0; row < 2; row++ ){ 
			for ( int col=0; col<2; col++ ){
				callCentral(col + row*nx);           // bottom edge
				callCentral(col + row* (nx+1) + nh); // left edge
			}

			for ( int col=2 ; col<nx-1; col++ ){
				callCentral(col + row*nx);         // bottom edge
				callMUSCL(col + row* (nx+1) + nh); // left edge
			}

			for ( int col=nx-1; col<nx; col++ ){
				callCentral(col + row*nx);           // bottom edge
				callCentral(col + row* (nx+1) + nh); // left edge
			}
			callCentral( nx + row*(nx+1) + nh ); // right edge
		}
		
		// loop through rows between first two rows and last two rows.
		OMP_PRAGMA(for schedule(static) nowait)
		for ( int row=2; row<ny-1; row++ ){
			for ( int col=0; col<2; col++ ){
				callMUSCL(col + row*nx);             // bottom edge
				callCentral(col + row* (nx+1) + nh); // left edge
			}

			for ( int col=2 ; col<nx-1; col++ ){
				callMUSCL(col + row*nx);           // bottom edge
				callMUSCL(col + row* (nx+1) + nh); // left edge
			}

			for ( int col=nx-1; col<nx; col++ ){
				callMUSCL(col + row*nx);             // bottom edge
				callCentral(col + row* (nx+1) + nh); // left edge
			}
			callCentral( nx + row*(nx+1) + nh ); // right edge
		}

		// last row
		
		OMP_PRAGMA(single nowait)
		{
			// first two cells
			for ( int col=0; col<2; col++ ){
				callCentral(col + (ny-1)*nx);          // bottom edge
				callCentral(col + (ny-1)*(nx+1) + nh); // left edge
				callCentral(col + ny*nx);              // top Edge
			}

			// middle cells
			for ( int col=2 ; col<nx-1; col++ ){
				callCentral(col + (ny-1)*nx);        // bottom edge
				callMUSCL(col + (ny-1)*(nx+1) + nh); // left edge
				callCentral(col + ny*nx);            // top Edge
				
			}
			//last cell
			callCentral(ny*nx - 1);          // bottom edge
			callCentral((ny+1)*nx -1);       // top edge
			callCentral(ny*(nx+1) + nh - 2); // left edge
			callCentral(ny*(nx+1) + nh - 1); // right edge
		}
	}
}

TRAVERSAL_FUNC_SECOND_ORDER(rowwiseTraversal_secondOrder_seamed){
	STRUCTURED_TRAVERSAL_ASSERT;
	auto callCentral = [&](Index e){
		funcCentral(mesh, e, std::forward<Ts>(args) ... );
	};
	auto callMUSCL = [&](Index e){
		funcMUSCL(mesh, e, std::forward<Ts>(args) ... );
	};

	Index
		nx { mesh.derived().size().x() },
		ny { mesh.derived().size().y() };
	
	Index nh { nHorizontalEdges( mesh.derived().size() ) };

	/** Info:
	 * threadsafety requires that adjacent rows must not be computed in parallel,
	 * since calculation on horizontal edges writes into the same array. 
	*/
	/* same procedure for the top row of cells and two bottom rows */
	OMP_PRAGMA(parallel)
	{
		OMP_PRAGMA(sections)
		{
			OMP_PRAGMA(section)
			{ // bottom rows, done by the same thread
				for ( int row = 0; row < 2; row++ ){
					for ( int col = 0; col < 2; col++ ){
						callCentral(col + row*nx); //bottom edge
						callCentral(col + row*(nx+1) + nh); // left edge
					}

					for ( int col = 2 ; col < nx-1; col++ ){
						callCentral(col + row*nx); // bottom edge
						callMUSCL(col + row*(nx+1) + nh); // left edge
					}

					//last cell
					callCentral((nx-1) + row*nx); //bottom edge
					callCentral((nx-1) + row*(nx+1) + nh); // left edge
					callCentral( nx + row*(nx+1) + nh );// right edge
				}
			}
			OMP_PRAGMA(section)
			{ // top row of cells done by one thread.
				for ( int col = 0; col < 2; col++ ){
					callCentral(col + (ny-1)*nx); //bottom edge
					callCentral(col + (ny-1)*(nx+1) + nh); // left edge
					callCentral(col + ny*nx); //top edge
				}

				for ( int col = 2 ; col < nx-1; col++ ){
					callCentral(col + (ny-1)*nx); // bottom edge
					callMUSCL(col + (ny-1)*(nx+1) + nh); // left edge
					callCentral(col + ny*nx); // top edge
				}

				// last cell
				callCentral( ny*nx - 1); //bottom edge
				callCentral(nh + ny*(nx+1) -2); // left edge
				callCentral(nh + ny*(nx+1) -1); // right edge
				callCentral( (ny+1)*nx -1); // top edge
			}
		} // implicit barrier
		/* work through the middle rows in parallel, 
		 * but always leave out every second row */
		OMP_PRAGMA(for schedule (static))
		for ( int row=2; row<ny-1; row+=2 ){
			for ( int col=0; col<2; col++ ){
				callMUSCL(col + row*nx); //bottom edge
				callCentral(col + row* (nx+1) + nh); // left edge
			}

			for ( int col=2 ; col<nx-1; col++ ){
				callMUSCL(col + row*nx); // bottom edge
				callMUSCL(col + row* (nx+1) + nh); // left edge
			}

			// last cell
			callMUSCL((nx-1) + row*nx); //bottom edge
			callCentral((nx-1) + row* (nx+1) + nh); // left edge
			callCentral( nx + row*(nx+1) + nh ); //right edge
		}
		/* work through rows that got left out */
		OMP_PRAGMA(for schedule (static) )
		for ( int row=3; row<ny-1; row+=2 ){
			for ( int col=0; col<2; col++ ){
				callMUSCL(col + row*nx); //bottom edge
				callCentral(col + row* (nx+1) + nh); // left edge
			}

			for ( int col=2 ; col<nx-1; col++ ){
				callMUSCL(col + row*nx); // bottom edge
				callMUSCL(col + row* (nx+1) + nh); // left edge
			}
			// last cell
			callMUSCL((nx-1) + row*nx); //bottom edge
			callCentral((nx-1) + row* (nx+1) + nh); // left edge
			callCentral( nx + row*(nx+1) + nh ); // right edge
		}
	}
}

/* A second order blockwise traversal that requires flux calculation 
 * functions that use the atomic update clause.
 */
TRAVERSAL_FUNC_SECOND_ORDER(blockwiseTraversal_secondOrder){
	STRUCTURED_TRAVERSAL_ASSERT;

	auto callCentral = [&](Index e){
		funcCentral(mesh, e, std::forward<Ts>(args) ... );
	};
	auto callMUSCL = [&](Index e){
		funcMUSCL(mesh, e, std::forward<Ts>(args) ... );
	};

	/* square root of number of cells per block */
	Index blocksize { mesh.derived().blocksize() };

	Index
		nx { mesh.derived().size().x() },
		ny { mesh.derived().size().y() },
		nh { nHorizontalEdges( mesh.derived().size() ) };

	Array2i nBlocks, rem, remBeg;
	for (int i=0; i<2; ++i){
		nBlocks[i] = mesh.derived().size()[i] / blocksize;
		rem    [i] = mesh.derived().size()[i] % blocksize;
	}
	
	// first cell of remainder section in x- and y-direction.
	remBeg[0] = nx - rem[0];
	remBeg[1] = ny - rem[1];

	if( rem.x() > 0 && rem.y() > 0){ 
		OMP_PRAGMA(parallel)
		{
			// inner blocks
			OMP_PRAGMA(for schedule(static) collapse(2) nowait)
			for(Index by = 1; by < nBlocks.y(); ++by) {
				for(Index bx = 1; bx < nBlocks.x(); ++bx) {
					for(Index row = by*blocksize; row < (by+1)*blocksize; ++row) {
						for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col) {
							callMUSCL(row*nx + col);
							callMUSCL(nh + row*(nx+1) + col);
						}
					}
				}
			}
			
			// left side
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index by = 1; by<nBlocks.y(); ++by){
				for(Index row =by*blocksize; row < (by+1)*blocksize; ++row){
					for(Index col = 0; col<2 ; ++col) {
						callMUSCL(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col<blocksize ; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
			}
			
			// bottom blocks
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index bx = 1; bx<nBlocks.x(); ++bx){
				for(Index row = 0; row < 2; ++row)
					for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col) {
						callCentral(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				for(Index row = 2; row < blocksize; ++row)
					for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
			}

			// right side
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index by = 1; by < nBlocks.y(); ++by) {
				for(Index row = by*blocksize; row < (by+1)*blocksize; ++row) {
					for(Index col = remBeg.x(); col < nx-1; ++col){
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
					callMUSCL((row+1)*nx - 1);
					callCentral( nh + (row+1)*(nx+1) -2 );
					callCentral( nh + (row+1)*(nx+1) -1 );
				}
			}
			
			// top side
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index bx = 1; bx < nBlocks.x(); ++bx){
				for(Index row = remBeg.y(); row < ny-1; ++ row){
					for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col){
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
				for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col){
						callCentral((ny-1)*nx + col);
						callCentral(ny*nx + col);
						callMUSCL(nh + (ny-1)*(nx+1) + col);
				}
			}

			
			// top left
			OMP_PRAGMA(single nowait)
			{
				for(Index row = remBeg.y(); row < ny-1; ++row){
					for(Index col = 0; col < 2; ++col){
						callMUSCL(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col < blocksize; ++col){
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
				// top row of block
				for(Index col = 0; col < 2; ++col){
					callCentral((ny-1)*nx + col);
					callCentral(ny*nx + col);
					callCentral(nh + (ny-1)*(nx+1) + col);

				}
				for(Index col = 2; col < blocksize; ++col){
					callCentral((ny-1)*nx + col);
					callCentral(ny*nx + col);
					callMUSCL(nh + (ny-1)*(nx+1) + col);
				}
			}
			// bottom right
			OMP_PRAGMA(single nowait)
			{
				for(Index row = 0; row<2; ++row){
					for(Index col = remBeg.x(); col < nx-1; ++col){
						callCentral(row * nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
					callCentral((row+1) * nx - 1 );
					callCentral(nh + (row+1)*(nx+1) -1);
					callCentral(nh + (row+1)*(nx+1) -2);
				}
				// rest of the rows have MUSCL horizontal edges
				for(Index row = 2; row<blocksize; ++row){
					for(Index col = remBeg.x(); col < nx-1; ++col){
						callMUSCL(row * nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
					callMUSCL  ((row+1) * nx - 1 );
					callCentral(nh + (row+1)*(nx+1) -2);
					callCentral(nh + (row+1)*(nx+1) -1);
				}
			}
			// top right corner block
			OMP_PRAGMA(single nowait)
			{
				for(Index row = remBeg.y(); row < ny-1; ++row){
					for(Index col = remBeg.x(); col < nx-1; ++col){
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
					callMUSCL((row+1)*nx - 1);
					callCentral(nh + (row+1)*(nx+1) - 2);
					callCentral(nh + (row+1)*(nx+1) - 1);
				}
				for(Index col = remBeg.x(); col< nx-1; ++col){
					callCentral((ny-1)*nx + col);
					callCentral(ny*nx + col);
					callMUSCL(nh + (ny-1)*(nx+1) + col);
				}
				callCentral(ny*nx-1);
				callCentral((ny+1)*nx-1);
				callCentral(nh + ny*(nx+1) - 2);
				callCentral(nh + ny*(nx+1) - 1);
			}

			// bottom left block
			OMP_PRAGMA(single nowait)
			{
				for(Index row = 0; row < 2; ++row){
					for(Index col = 0; col < 2; ++col){
						callCentral(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col < blocksize; ++col) {
						callCentral(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
				for(Index row = 2; row < blocksize; ++row){
					for(Index col = 0; col < 2; ++col){
						callMUSCL(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col < blocksize; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
			}

		}
	}

	// only a right remainder
	else if (rem.x() > 0 && rem.y() == 0){
		OMP_PRAGMA(parallel)
		{
			// inner real blocks
			OMP_PRAGMA(for schedule(static) collapse(2) nowait)
			for(Index by = 1; by < nBlocks.y()-1; ++by) {
				for(Index bx = 1; bx < nBlocks.x(); ++bx) {
					for(Index row = by*blocksize; row < (by+1)*blocksize; ++row) {
						for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col) {
							callMUSCL(row*nx + col);
							callMUSCL(nh + row*(nx+1) + col);
						}
					}
				}
			}
			
			// left side
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index by = 1; by<nBlocks.y()-1; ++by){
				for(Index row =by*blocksize; row < (by+1)*blocksize; ++row){
					for(Index col = 0; col<2 ; ++col) {
						callMUSCL(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col<blocksize ; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
			}
			// bottom blocks
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index bx = 1; bx<nBlocks.x(); ++bx){
				for(Index row = 0; row < 2; ++row)
					for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col) {
						callCentral(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				for(Index row = 2; row < blocksize; ++row)
					for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
			}

			// right side
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index by = 1; by < nBlocks.y()-1; ++by) {
				for(Index row = by*blocksize; row < (by+1)*blocksize; ++row) {
					for(Index col = remBeg.x(); col < nx-1; ++col){
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
					callMUSCL((row+1)*nx - 1);
					callCentral( nh + (row+1)*(nx+1) -2 );
					callCentral( nh + (row+1)*(nx+1) -1 );
				}
			}

			// top side
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index bx = 1; bx<nBlocks.x(); ++bx){
				for(Index row = ny-blocksize; row < ny-1; ++row){
					for(Index col = bx*blocksize; col<(bx+1)*blocksize; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
				for(Index col = bx*blocksize; col<(bx+1)*blocksize; ++col) {
					callCentral((ny-1)*nx + col);
					callCentral(ny*nx + col);
					callMUSCL(nh + (ny-1)*(nx+1) + col);
				}
			}

			// bottom right			
			OMP_PRAGMA(single nowait)
			{
				for(Index row = 0; row<2; ++row){
					for(Index col = remBeg.x(); col < nx-1; ++col){
						callCentral(row * nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
					callCentral((row+1) * nx - 1 );
					callCentral(nh + (row+1)*(nx+1) -1);
					callCentral(nh + (row+1)*(nx+1) -2);
				}
				// rest of the rows have MUSCL horizontal edges
				for(Index row = 2; row<blocksize; ++row){
					for(Index col = remBeg.x(); col < nx-1; ++col){
						callMUSCL(row * nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
					callMUSCL  ((row+1) * nx - 1 );
					callCentral(nh + (row+1)*(nx+1) -2);
					callCentral(nh + (row+1)*(nx+1) -1);
				}
			}
			
			// top right 
			OMP_PRAGMA(single nowait)
			{
				for(Index row = ny-blocksize; row < ny-1; ++row){
					for(Index col = remBeg.x(); col < nx-1; ++col) {
						callMUSCL(row * nx + col);
						callMUSCL(nh + row * (nx+1) + col);
					}
					callMUSCL((row+1)*nx - 1);
					callCentral(nh + (row+1)*(nx+1) - 2);
					callCentral(nh + (row+1)*(nx+1) - 1);
				}
				for(Index col = remBeg.x(); col < nx-1; ++col) {
					callCentral((ny-1) * nx + col);
					callCentral(ny*nx + col);
					callMUSCL(nh + (ny-1) * (nx+1) + col);
				}
				callCentral((ny+1)*nx - 1);
				callCentral(ny*nx - 1);
				callCentral(nh + (ny)*(nx+1) - 2);
				callCentral(nh + (ny)*(nx+1) - 1);
			}

			// bottom left
			OMP_PRAGMA(single nowait)
			{
				for(Index row = 0; row < 2; ++row){
					for(Index col = 0; col < 2; ++col){
						callCentral(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col < blocksize; ++col) {
						callCentral(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
				for(Index row = 2; row < blocksize; ++row){
					for(Index col = 0; col < 2; ++col){
						callMUSCL(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col < blocksize; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
			}

			// top left 
			OMP_PRAGMA(single nowait)
			{
				for(Index row = ny-blocksize; row < ny-1; ++row){
					for(Index col = 0; col<2; ++col) {
						callMUSCL(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col<blocksize; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
				for(Index col = 0; col<2; ++col) {
					callCentral((ny-1)*nx + col);
					callCentral(ny*nx + col);
					callCentral(nh + (ny-1)*(nx+1) + col);
				}
				for(Index col = 2; col<blocksize; ++col) {
					callCentral((ny-1)*nx + col);
					callCentral(ny*nx + col);
					callMUSCL(nh + (ny-1)*(nx+1) + col);
				}
			}
			}

	}

	// only remainder on top
	else if(rem.x() == 0 && rem.y() > 0){
		OMP_PRAGMA(parallel)
		{
			// inner blocks
			OMP_PRAGMA(for schedule(static) collapse(2) nowait)
			for(Index by = 1; by < nBlocks.y(); ++by) {
				for(Index bx = 1; bx < nBlocks.x()-1; ++bx) {
					for(Index row = by*blocksize; row < (by+1)*blocksize; ++row) {
						for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col) {
							callMUSCL(row*nx + col);
							callMUSCL(nh + row*(nx+1) + col);
						}
					}
				}
			}
		
			// left side
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index by = 1; by<nBlocks.y(); ++by){
				for(Index row =by*blocksize; row < (by+1)*blocksize; ++row){
					for(Index col = 0; col<2 ; ++col) {
						callMUSCL(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col<blocksize ; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
			}

			// bottom blocks
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index bx = 1; bx<nBlocks.x()-1; ++bx){
				for(Index row = 0; row < 2; ++row)
					for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col) {
						callCentral(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				for(Index row = 2; row < blocksize; ++row)
					for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
			}
			
			// blocks on right side
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index by = 1; by<nBlocks.y(); ++by){
				for(Index row = by*blocksize; row < (by+1)*blocksize; ++row){
					for(Index col = nx-blocksize; col < nx-1; ++col){
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
					callMUSCL((row+1)*nx - 1);
					callCentral(nh + (row+1)*(nx+1) - 2);
					callCentral(nh + (row+1)*(nx+1) - 1);
				}
			}

			// top blocks
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index bx = 1; bx < nBlocks.x()-1; ++bx){
				for(Index row = remBeg.y(); row < ny-1; ++ row){
					for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col){
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
				for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col){
						callCentral((ny-1)*nx + col);
						callCentral(ny*nx + col);
						callMUSCL(nh + (ny-1)*(nx+1) + col);
				}
			}

			// bottom right block
			OMP_PRAGMA(single nowait)
			{
				for(Index row = 0; row < 2; ++row){
					for(Index col = nx-blocksize; col< nx-1; ++col) {
						callCentral(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
					callCentral((row+1)*nx - 1);
					callCentral(nh + (row+1)*(nx+1) - 2);
					callCentral(nh + (row+1)*(nx+1) - 1);
				}
				for(Index row = 2; row < blocksize; ++row){
					for(Index col = nx-blocksize; col< nx-1; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
					callMUSCL((row+1)*nx - 1);
					callCentral(nh + (row+1)*(nx+1) - 2);
					callCentral(nh + (row+1)*(nx+1) - 1);
				}
			}
			
			// top left
			OMP_PRAGMA(single nowait)
			{
				for(Index row = remBeg.y(); row < ny-1; ++row){
					for(Index col = 0; col < 2; ++col){
						callMUSCL(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col < blocksize; ++col){
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
				// top row of block
				for(Index col = 0; col < 2; ++col){
					callCentral((ny-1)*nx + col);
					callCentral(ny*nx + col);
					callCentral(nh + (ny-1)*(nx+1) + col);

				}
				for(Index col = 2; col < blocksize; ++col){
					callCentral((ny-1)*nx + col);
					callCentral(ny*nx + col);
					callMUSCL(nh + (ny-1)*(nx+1) + col);
				}
			}
			
			// top right
			OMP_PRAGMA(single nowait)
			{
				for(Index row = remBeg.y(); row < ny-1; ++row){
					for(Index col = nx-blocksize; col <nx-1; ++col){
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
					callMUSCL((row+1)*nx -1);
					callCentral(nh + (row+1)*(nx+1) - 2);
					callCentral(nh + (row+1)*(nx+1) - 1);
				}
				for(Index col = nx-blocksize; col <nx-1; ++col){
					callCentral((ny-1)*nx + col);
					callCentral(ny*nx + col);
					callMUSCL(nh + (ny-1)*(nx+1) + col);
				}
				callCentral(ny*nx -1);
				callCentral((ny+1)*nx - 1);
				callCentral(nh + ny*(nx+1) - 2);
				callCentral(nh + ny*(nx+1) - 1);
			}
			
			// bottom left
			OMP_PRAGMA(single nowait)
			{
				for(Index row = 0; row < 2; ++row){
					for(Index col = 0; col < 2; ++col){
						callCentral(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col < blocksize; ++col) {
						callCentral(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
				for(Index row = 2; row < blocksize; ++row){
					for(Index col = 0; col < 2; ++col){
						callMUSCL(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col < blocksize; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
			}

			}
	}

	// no remainder block
	else {
		OMP_PRAGMA(parallel)
		{

			// inner blocks
			OMP_PRAGMA(for schedule(static) collapse(2) nowait)
			for(Index by = 1; by < nBlocks.y()-1; ++by) {
				for(Index bx = 1; bx < nBlocks.x()-1; ++bx) {
					for(Index row = by*blocksize; row < (by+1)*blocksize; ++row) {
						for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col) {
							callMUSCL(row*nx + col);
							callMUSCL(nh + row*(nx+1) + col);
						}
					}
				}
			}

			// left blocks
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index by = 1; by<nBlocks.y()-1; ++by){
				for(Index row =by*blocksize; row < (by+1)*blocksize; ++row){
					for(Index col = 0; col<2 ; ++col) {
						callMUSCL(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col<blocksize ; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
			}

			// bottom blocks
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index bx = 1; bx<nBlocks.x()-1; ++bx){
				for(Index row = 0; row < 2; ++row)
					for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col) {
						callCentral(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				for(Index row = 2; row < blocksize; ++row)
					for(Index col = bx*blocksize; col < (bx+1)*blocksize; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
			}

			// right blocks
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index by = 1; by<nBlocks.y()-1; ++by){
				for(Index row = by*blocksize; row < (by+1)*blocksize; ++row){
					for(Index col = nx-blocksize; col < nx-1; ++col){
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
					callMUSCL((row+1)*nx - 1);
					callCentral(nh + (row+1)*(nx+1) - 2);
					callCentral(nh + (row+1)*(nx+1) - 1);
				}
			}
			
			// top blocks
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index bx = 1; bx<nBlocks.x()-1; ++bx){
				for(Index row = ny-blocksize; row < ny-1; ++row){
					for(Index col = bx*blocksize; col<(bx+1)*blocksize; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
				for(Index col = bx*blocksize; col<(bx+1)*blocksize; ++col) {
					callCentral((ny-1)*nx + col);
					callCentral(ny*nx + col);
					callMUSCL(nh + (ny-1)*(nx+1) + col);
				}
			}
			
			// bottom right
			OMP_PRAGMA(single nowait)
			{
				for(Index row = 0; row < 2; ++row){
					for(Index col = nx-blocksize; col< nx-1; ++col) {
						callCentral(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
					callCentral((row+1)*nx - 1);
					callCentral(nh + (row+1)*(nx+1) - 2);
					callCentral(nh + (row+1)*(nx+1) - 1);
				}
				for(Index row = 2; row < blocksize; ++row){
					for(Index col = nx-blocksize; col< nx-1; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
					callMUSCL((row+1)*nx - 1);
					callCentral(nh + (row+1)*(nx+1) - 2);
					callCentral(nh + (row+1)*(nx+1) - 1);
				}
			}
			// top left
			OMP_PRAGMA(single nowait)
			{
				for(Index row = ny-blocksize; row < ny-1; ++row){
					for(Index col = 0; col<2; ++col) {
						callMUSCL(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col<blocksize; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
				for(Index col = 0; col<2; ++col) {
					callCentral((ny-1)*nx + col);
					callCentral(ny*nx + col);
					callCentral(nh + (ny-1)*(nx+1) + col);
				}
				for(Index col = 2; col<blocksize; ++col) {
					callCentral((ny-1)*nx + col);
					callCentral(ny*nx + col);
					callMUSCL(nh + (ny-1)*(nx+1) + col);
				}
			}
			// top right
			OMP_PRAGMA(single nowait)
			{
				for(Index row = ny-blocksize; row < ny-1; ++row){
					for(Index col = nx-blocksize; col < nx-1; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
					callMUSCL((row+1)*nx - 1);
					callCentral(nh + (row+1)*(nx+1) - 2);
					callCentral(nh + (row+1)*(nx+1) - 1);
				}
				for(Index col = nx-blocksize; col < nx-1; ++col){
					callCentral((ny-1)*nx + col);
					callCentral( ny   *nx + col);
					callMUSCL(nh + (ny-1)*(nx+1) + col);
				}
				callCentral(ny*nx - 1);
				callCentral((ny+1)*nx - 1);
				callCentral(nh + ny*(nx+1) - 2);
				callCentral(nh + ny*(nx+1) - 1);
			}
			// bottom left 
			OMP_PRAGMA(single nowait)
			{
				for(Index row = 0; row < 2; ++row){
					for(Index col = 0; col < 2; ++col){
						callCentral(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col < blocksize; ++col) {
						callCentral(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
				for(Index row = 2; row < blocksize; ++row){
					for(Index col = 0; col < 2; ++col){
						callMUSCL(row*nx + col);
						callCentral(nh + row*(nx+1) + col);
					}
					for(Index col = 2; col < blocksize; ++col) {
						callMUSCL(row*nx + col);
						callMUSCL(nh + row*(nx+1) + col);
					}
				}
			}

		}
	}
}


TRAVERSAL_FUNC_SECOND_ORDER(blockwiseTraversal_secondOrder_seamed){
	STRUCTURED_TRAVERSAL_ASSERT;

	auto callCentral = [&](Index e){
		funcCentral(mesh, e, std::forward<Ts>(args) ... );
	};
	auto callMUSCL = [&](Index e){
		funcMUSCL(mesh, e, std::forward<Ts>(args) ... );
	};

	Index blocksize { mesh.derived().blocksize() };

	Index
		nx { mesh.derived().size().x() },
		ny { mesh.derived().size().y() },
		nh { nHorizontalEdges( mesh.derived().size() ) },
		nEdges { mesh.nEdges() };
	
	Array2i innerLength { nx-2, ny-2 };
	Array2i nBlocks, rem, remBeg;
	for(int i = 0; i<2; i++){
		nBlocks[i] = innerLength[i] / blocksize;
		rem[i] = innerLength[i] % blocksize;
	}
	remBeg[0] = nx - rem[0] - 1;
	remBeg[1] = ny - rem[1] - 1;


	OMP_PRAGMA(parallel )
	{

		/* inner blocks */
		OMP_PRAGMA(for schedule(static) collapse(2) nowait)
		for(Index by = 0; by<nBlocks.y(); ++by){
			for(Index bx = 0; bx<nBlocks.x(); ++bx){
				for(Index row = 0; row<blocksize-1; row++){
					
					Index begHori { nx+1 +    by*blocksize*nx +     row*nx +     bx*blocksize };
					Index begVert { nh+nx+2 + by*blocksize*(nx+1) + row*(nx+1) + bx*blocksize };

					callMUSCL(begHori + nx);
					for(Index col = 1; col<blocksize; col++){
						callMUSCL(begHori + nx + col);
						callMUSCL(begVert + col);
					}
				}
				for(Index col = 1; col<blocksize; col++){
					callMUSCL(nh+nx+2 + (((by+1)*blocksize)-1)*(nx+1) + bx*blocksize + col );
				}
			}
		}
		// remainder right
		if(rem.x() >= 1){
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for(Index by = 0; by<nBlocks.y(); by++){ // skip last block
				
				for(Index row = 0 ; row<blocksize-1; row++){
				
					// first left and bottom edge of current row in remainder block.
					Index begHori { nx +      by*blocksize*nx +     row*nx     + remBeg.x() };
					Index begVert { nh+nx+1 + by*blocksize*(nx+1) + row*(nx+1) + remBeg.x() };

					callMUSCL( begHori + nx );

					for(Index col = 1; col < rem.x(); col++){
						
						callMUSCL( begHori + nx + col ); // horizontal on top
						callMUSCL( begVert + col ); // vertical to the left.
					}
				}

				// top row of a remainder block only has vertical edges
				for(Index col = 1; col < rem.x(); col++){
					callMUSCL( (nh+nx+1) + ((by+1)*blocksize-1)*(nx+1) + remBeg.x() + col );
				}
			}
		}

		// remainder top
		if(rem.y() >= 1){
			OMP_PRAGMA(for schedule(dynamic) nowait)
			for( Index bx = 0; bx<nBlocks.x(); bx++){
				
				Index begVert { nh + remBeg.y()*(nx+1)+1 + bx*blocksize };

				for(Index col = 1; col < blocksize; col++){
					callMUSCL(begVert + col);
				}

				for(Index row = 1; row<rem.y(); row++){

					// first left and bottom edge of current row in remainder block.
					Index begHori { remBeg.y()*nx+1 +          row*nx +     bx*blocksize };
					begVert =       nh + remBeg.y()*(nx+1)+1 + row*(nx+1) + bx*blocksize;
					
					callMUSCL( begHori); // first horizontal, no vertical to the left.
					for(Index col = 1; col<blocksize; col++){
						callMUSCL( begHori + col);
						callMUSCL( begVert + col);
					}
				}
			}
		}

		// remainder top right corner
		OMP_PRAGMA(single)
		{
			if(rem.x() > 1 && rem.y() > 1){

				Index begVert { nh + remBeg.y()*(nx+1) + remBeg.x() };
				Index begHori { remBeg.y()*nx + remBeg.x() };

				for(Index col = 1; col<rem.x(); ++col) {
					callMUSCL(begVert+col);
				}

				for ( Index row = 1; row < rem.y(); ++row){
					callMUSCL( begHori + nx*row );

					for ( Index col = 1; col< rem.x(); ++col){
						callMUSCL( begHori + row*nx + col);
						callMUSCL( begVert + row*(nx+1)+ col );
					}
				}
			}

			else if(rem.x() > 1 && rem.y() == 1) {

				Index begVert { nh + remBeg.y()*(nx+1) + remBeg.x() };

				for ( Index col = 1; col< rem.x(); ++col)
					callMUSCL (begVert+col);
			}

			else if(rem.x() == 1 && rem.y() > 1){
				
				Index begHori { remBeg.y()*nx + remBeg.x() };
				
				for( Index row = 1; row < rem.y() ; ++row)
					callMUSCL( begHori + row*nx );
			}
		} // implicit barrier

		// vertical Central edges in left and right cell columns
		// start with lower left in a single
		OMP_PRAGMA(single nowait)
		{
			callCentral(nh);
			callCentral(nh+1);
		}
		
		// loop treats last two edges of a row and first two edges of second row
		// in hope of faster memory access since these 4 edges are adjacently stored.
		OMP_PRAGMA(for schedule(static) nowait)
		for(Index row = 1; row<ny; row++){ 
			callCentral(nh + row*(nx+1)-2);
			callCentral(nh + row*(nx+1)-1);
			callCentral(nh + row*(nx+1));
			callCentral(nh + row*(nx+1)+1);
		}
		
		// last two vertical edges, top right.
		OMP_PRAGMA(single nowait)
		{
			callCentral(nEdges-1);
			callCentral(nEdges-2);
		}

		// shared vertical MUSCL edges
		OMP_PRAGMA(for schedule(static) collapse(2))
		for(Index col = 1; col < ((rem.x() > 0) ? nBlocks.x()+1 : nBlocks.x() ) ; col++){
			for(Index row = 0; row<innerLength.y(); row++){
				callMUSCL(nh+nx+2 + row*(nx+1) + col*blocksize );
			}
		}

		/* MUSCL scheme for horizontal edges in left and right column */
		// first horizontal edge
		OMP_PRAGMA(barrier)

		OMP_PRAGMA(single)
		{
			callMUSCL(2*nx);
		}
		// every other
		OMP_PRAGMA(for schedule(static))
		for(Index row=3; row<ny-1; row += 2){
			callMUSCL(row*nx-1);
			callMUSCL(row*nx);
		}
		// every skipped 
		OMP_PRAGMA(for schedule(static))
		for(Index row=4; row<ny-1; row += 2){
			callMUSCL(row*nx-1);
			callMUSCL(row*nx);
		}
		// last horizontal top-right
		OMP_PRAGMA(single )
		{
			callMUSCL((ny-1)*nx - 1);
		}

		/* central scheme on horizontal edges */
		OMP_PRAGMA(for schedule(static) nowait)
		for(Index col = 0; col<nx; col++){
			callCentral(col);
			callCentral(nx+col);
			callCentral((ny-1)*nx+col);
			callCentral(ny*nx+col);
		}
		
		// shared horizontal MUSCL edges
		OMP_PRAGMA(for schedule(static) collapse(2))
		for(Index row = 1; row < ((rem.y() > 0) ? nBlocks.y()+1 : nBlocks.y() ) ; row++){
			for(Index col = 0; col<innerLength.x(); col++){
				callMUSCL(nx+1 + row*blocksize*nx + col);
			}
		}
		/* MUSCL scheme for vertical edges in bottom and top row*/
		// every second edge, for thread safety
		// cache friendly for larger meshes.
		OMP_PRAGMA(for schedule(static))
		for(Index col = 2; col<nx-1; col += 2){
			// Again, maybe split the loop into two loops for better memory access.
			callMUSCL(nh + col); // bottom row
			callMUSCL(nEdges - (nx+1) + col); // top row
		}

		// every other edge
		OMP_PRAGMA(for schedule(static))
		for(Index col = 3; col<nx-1; col += 2){
			// Again, maybe split the loop into two loops for better memory access.
			callMUSCL(nh + col); // bottom row
			callMUSCL(nEdges - (nx+1) + col); // top row
		}
	}
}

#undef FORWARD_ARGS

} // end namespace hms